bulk-org
Raphael Kubler
7/26/2023
Introduction
Bulk RNA-seq analysis of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).
CO vs COiMG
Organoids co-cultured with vs without microglia.
Variance analysis
PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
theme_classic()
Heatmap
#for control only
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'
#horizontal
ra <- rowAnnotation(
COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
col = list(
COiMG = c("COiMg"='tomato','CO'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr <- t(batch.rem3_mod1[,meta3$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
DEA
Volcano
volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
res <-
mutate(res,
sig = case_when(
padj >= 0.05 ~ "non_sig",
padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
)) %>%
mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
mutate(log2FoldChange = case_when(
log2FoldChange > 10 ~ Inf,
log2FoldChange < -10 ~ -Inf,
TRUE ~ log2FoldChange
)) %>%
mutate(class = paste(sig, direction))
if( type == "ALS"){
xpos <- 0.5
ymax <- ymax1
xlim <- c(xmax1,xmax2)
}else{
xpos <- 0.025
ymax <- ymax2
xlim <- c(-0.042, 0.042)
}
de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
filter(sig != "non_sig") %>%
mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
mutate(n = formatC(n, format="f", big.mark=",", digits=0))
plot <- res %>%
mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) +
#geom_point(aes(colour = class ), size = 0.5) +
rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
scale_colour_manual(values = c("non_sig up" = "gray",
"non_sig down" = "gray",
"sig up" = "#EB7F56",
"sig - strong up" = "#B61927",
"sig down" = "#4F8FC4",
"sig - strong down" = "dodgerblue4"
)) +
theme_bw() +
labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
guides(colour = FALSE) +
scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.border = element_blank(),
axis.ticks = element_line(colour = "black")
) +
geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
scale_x_continuous(limits = xlim)
if(!is.null(annotate_by)){
plot <- plot +
ggrepel::geom_text_repel(
fontface = "italic",
data = filter(res, symbol %in% annotate_by),
aes(x = log2FoldChange, y = -log10(pvalue), label = symbol),
min.segment.length = unit(0, "lines"),
size = 2.3) +
geom_point(
data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
) +
geom_point(aes(colour = class ),
data = filter(res, symbol %in% annotate_by), size = 0.6
)
}
return(plot)
}
res_COiMG$symbol <- rownames(res_COiMG)
volcano_plot(data.frame(res_COiMG), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top50.coimg <- rbind(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names <- rownames(top50.coimg)
ra2 <- rowAnnotation(
COiMG = meta3$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra2,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
coimg.deg <- rbind(results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,1]),],
results$COiMG[rownames(results$COiMG) %in% na.omit(all_res$COiMG$DEGs[,2]),])
df <- cbind('DEG'=rownames(coimg.deg),'lFC'=round(coimg.deg$log2FoldChange,2),'Cilia-associated'=rownames(coimg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)GO
dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$COiMG$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts3)]), check.names = F)
for(i in names(cgenes)[-1]){
k <- na.omit(cgenes[i])
s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts3),]
m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
IFNy stimulation
Let’s do this for IFNy stimulation
Variance analysis
PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
#stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
theme_classic()
DEA
Volcano
res_IFNg$symbol <- rownames(res_IFNg)
volcano_plot(data.frame(res_IFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
COiMg.degs <- all_res$COiMG$DEGsTop 50 DEGs heatmap
top50.ifng <- rbind(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names2 <- rownames(top50.ifng)
ra3 <- rowAnnotation(
Stimulation = meta3$condition,
col = list(
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra3,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
IFNg.deg <- rbind(results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,1]),],
results$IFNg[rownames(results$IFNg) %in% na.omit(all_res$IFNg$DEGs[,2]),])
df <- cbind('DEG'=rownames(IFNg.deg),'lFC'=round(IFNg.deg$log2FoldChange,2),'Cilia-associated'=rownames(IFNg.deg) %in% organoid$CBC)
createDT(data.frame(df), "Status", scrollY=1000)GO
dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
For organoid cell-types.
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For NDD gene lists.
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]## [1] "CYP7B1" "GRIN2A"For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
Cell-type marker analysis
Cell-types with only 1 gene found back in our expression data carry principal component expression = 0.
Heatmaps of gene & PC expression
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
cell.types <- as.matrix(readxl::read_xlsx('Typical markers cell types organoids 07182023.xlsx', col_names = F))
PC.covs <- matrix(ncol=length(cell.types[,1]),nrow=length(colnames(batch.rem3_mod1)))
colnames(PC.covs) <- cell.types[,1]
for(i in cell.types[,1]){
pca.cov <- prcomp(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],]))
PC.covs[,i] <- pca.cov$x[,1]
#if (table(rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1])[[2]] == 1){
# PC.covs[,i] <- t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],])
#}
}
rownames(PC.covs) <- colnames(batch.rem3_mod1)
ra4 <- rowAnnotation(
Stimulation = meta3$condition,
COiMG = meta3$COiMg,
col = list(COiMG = c('yes'='black','no'='white'),
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(PC.covs,
cluster_rows = F,
cluster_columns = F,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra4,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(PC.covs), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
#Let's do that for ctr/ifng and coimg/CO only
#for control only
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
#horizontal
ctr.ra <- rowAnnotation(
COiMG = ctr.meta$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr.pca <- PC.covs[meta3$condition %in% 'ctrl',]
ComplexHeatmap::Heatmap(ctr.pca,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ctr.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
ifng.meta <- meta3[!meta3$condition %in% 'LPS',]
#horizontal
ifng.ra <- rowAnnotation(
COiMG = ifng.meta$COiMg,
Stimulation = ifng.meta$condition,
col = list(COiMG = c('yes'='black','no'='white'),
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ifng.pca <- PC.covs[!meta3$condition %in% 'LPS',]
ComplexHeatmap::Heatmap(ifng.pca,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ifng.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
#Let's check only choroid plexus & ependymal
cp.ep <- c(na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][1,][-1]),
na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][2,][-1]))
ctr.df <- batch.rem3_mod1[,colnames(batch.rem3_mod1) %in% rownames(ctr.meta)]
ComplexHeatmap::Heatmap(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ctr.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,]))), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Boxplots
Let’s make boxplots instead
ifng.pca <- PC.covs[!meta3$condition %in% 'LPS',]
ifng.pca_long <- reshape::melt(ifng.pca)
colnames(ifng.pca_long) <- c('Sample','Cell.type','PC')
ifng.pca_long$Stimulation <- ifng.meta$condition
ggplot(ifng.pca_long,aes_string(x='Cell.type',y='PC', fill='Stimulation')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
ctr.pca <- PC.covs[meta3$condition %in% 'ctrl',]
coimg.pca_long <- reshape::melt(ctr.pca)
colnames(coimg.pca_long) <- c('Sample','Cell.type','PC')
coimg.pca_long$Coculture <- ctr.meta$COiMg
ggplot(coimg.pca_long,aes_string(x='Cell.type',y='PC', fill='Coculture')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
Dissecting DEA results
CO-IFNy
These are the results of the DEA of CTR vs IFNy in CO (not co-cultured with microglia) only. Cilia-associated genes remain downregulated. We find this as well in CO vs COiMG and more so also in our single-cell data, which did not have any IFNy-stimulated samples. Meaning: this is a result robust to IFNy AND microglia co-culture separately.
load('/Users/kubler01/Documents/R projects/bulk-org/07262023_IFNy in CO.RData')
new_IFNy$padj[is.na(new_IFNy$padj)] <- 1
downreg.genes <- rownames(new_IFNy[new_IFNy$log2FoldChange < -1 & new_IFNy$padj < 0.05,])
upreg.genes <- rownames(new_IFNy[new_IFNy$log2FoldChange > 1 & new_IFNy$padj < 0.05,])
IFNy.res <- new_IFNy[rownames(new_IFNy) %in% c(upreg.genes,downreg.genes),]
df <- cbind('DEG'=rownames(IFNy.res),'lFC'=round(IFNy.res$log2FoldChange,2),'Cilia-associated'=rownames(IFNy.res) %in% organoid$CBC)
createDT(data.frame(df), "Stimulation", scrollY=1000)Volcano
new_IFNy$symbol <- rownames(new_IFNy)
volcano_plot(data.frame(new_IFNy), title = 'ctrl vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
CO-COiMg
This is to show that the DEA effect of cilia-associated genes remains present in our bulk-seq analysis when we subset to control samples only and test CO vs COiMg.
setwd('/Users/kubler01/Documents/R projects/bulk-org')
load('07262023_CO vs COiMG.RData')
new_coimg$padj[is.na(new_coimg$padj)] <- 1
downreg.genes <- rownames(new_coimg[new_coimg$log2FoldChange < -1 & new_coimg$padj < 0.05,])
upreg.genes <- rownames(new_coimg[new_coimg$log2FoldChange > 1 & new_coimg$padj < 0.05,])
coimg.res <- new_coimg[rownames(new_coimg) %in% c(upreg.genes,downreg.genes),]
df <- cbind('DEG'=rownames(coimg.res),'lFC'=round(coimg.res$log2FoldChange,2),'Cilia-associated'=rownames(coimg.res) %in% organoid$CBC)
createDT(data.frame(df), "COiMG", scrollY=1000)Volcano
new_coimg$symbol <- rownames(new_coimg)
volcano_plot(data.frame(new_coimg), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Correlation of lFCs
overlap <- rownames(coimg.res[rownames(coimg.res) %in% rownames(IFNy.res),])
cor(coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange)## [1] 0.8245697createDT(cbind('Genes'=overlap,
'COiMg'=coimg.res[rownames(coimg.res) %in% overlap,]$log2FoldChange,
'IFNy'=IFNy.res[rownames(IFNy.res) %in% overlap,]$log2FoldChange), "COiMG", scrollY=1000)